Anne Badel et Frédéric Guyon
2019-02-19
Y a-t-il des groupes ?
Oui, 4 groupes.
Trois grands principes de méthodes basées sur:
En fait, trois façons de voir les mêmes algorithmes
On considère les données comme des points de \(R^n\):
Ces données sont un classique des méthodes d'apprentissage
Dans un premier temps, regardons les données
dim(mes.iris)[1] 150 4
head(mes.iris) Sepal.Length Sepal.Width Petal.Length Petal.Width
1 5.1 3.5 1.4 0.2
2 4.9 3.0 1.4 0.2
3 4.7 3.2 1.3 0.2
4 4.6 3.1 1.5 0.2
5 5.0 3.6 1.4 0.2
6 5.4 3.9 1.7 0.4
On peut ensuite essayer de visualiser les données
plot(mes.iris)Sur la base d'une distance (souvent euclidienne)
Définition d'une distance : fonction positive de deux variables
Si 1,2,3 : dissimilarité
distance euclidienne ou distance \(l_2\): \(d(x,y)=\sqrt{\sum_i (x_i-y_i)^2}\)
distance de manahattan ou distance \(l_1\): \(d(x,y)=\sum_i |x_i-y_i|\)
distance du maximum ou l-infini, \(l_\infty\): \(d(x,y)=\max_i |x_i-y_i|\)
distance de Minkowski \(l_p\): \[d(x,y)=\sqrt[p]{\sum_i (|x_i-y_i|^p}\]
distance de Canberra (x et y valeurs positives): \[d(x,y)=\sum_i \frac{x_i-y_i}{x_i+y_i}\]
distance binaire ou distance de Jaccard ou Tanimoto: proportion de propriétés communes
Utilisées en bio-informatique:
\[d("BONJOUR", "BONSOIR")=2\]
Comme vu lors de la séance 3, il y a d'autres mesures de distances :
ne sont pas des distances, mais indices de dissimilarité:
rq : lors du TP, sur les données d'expression RNA-seq, nous utiliserons le coefficient de corrélation de Spearman et la distance associée, \(d = 1-r^2\)
\(D(C_1,C_2) = \min_{i \in C_1, j \in C_2} D(x_i, x_j)\)
\(D(C_1,C_2) = \max_{i \in C_1, j \in C_2} D(x_i, x_j)\)
\(D(C_1,C_2) = \frac{1}{N_1 N_2} \sum_{i \in C_1, j \in C_2} D(x_i, x_j)\)
\(d^2(C_i,C_j)= I_{intra}(C_i \cup C_j)-I_{intra}(C_i)-I_{intra}(C_j)\) \(D(C_1,C_2) = \sqrt{\frac{N_1N_2}{N_1 + N_2}} \| m_1 -m_2 \|\)
Ces données sont un classique des méthodes d'apprentissage
Dans un premier temps, regardons les données
dim(mes.iris)[1] 150 4
head(mes.iris) Sepal.Length Sepal.Width Petal.Length Petal.Width
1 5.1 3.5 1.4 0.2
2 4.9 3.0 1.4 0.2
3 4.7 3.2 1.3 0.2
4 4.6 3.1 1.5 0.2
5 5.0 3.6 1.4 0.2
6 5.4 3.9 1.7 0.4
str(mes.iris)'data.frame': 150 obs. of 4 variables:
$ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
$ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
$ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
$ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
summary(mes.iris) Sepal.Length Sepal.Width Petal.Length Petal.Width
Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
Median :5.800 Median :3.000 Median :4.350 Median :1.300
Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
On peut ensuite essayer de visualiser les données
plotplot(mes.iris)imageimage(1:nb.var, 1:nb.iris ,t(as.matrix(mes.iris)))Avant de commencer à travailler, il est nécessaire de commencer par vérifier que :
sum(is.na(mes.iris))[1] 0
iris.var <- apply(mes.iris, 2, var)
kable(iris.var, digits = 3)| x | |
|---|---|
| Sepal.Length | 0.686 |
| Sepal.Width | 0.190 |
| Petal.Length | 3.116 |
| Petal.Width | 0.581 |
sum(apply(mes.iris, 2, var) == 0)[1] 0
Afin de pouvoir considérer que toutes les variables sont à la même échelle, il est parfois nécessaire de normaliser les données.
mes.iris.centre <- scale(mes.iris, center=TRUE, scale=FALSE)mes.iris.scaled <- scale(mes.iris, center=TRUE, scale=TRUE)par(mfrow=c(1,2))
plot(mes.iris)plot(mes.iris.scaled)
par(mfrow=c(1,1)) ! ne pas faire si "grosses" données
par(mfrow=c(1,2))
image(1:nb.var, 1:nb.iris, t(as.matrix(mes.iris)), main="données brutes")
image(1:nb.var, 1:nb.iris, t(as.matrix(mes.iris.scaled)), main="données normalisées")par(mfrow=c(1,1))par(mfrow=c(1,2))
biplot(prcomp(mes.iris), main="données non normalisées")
biplot(prcomp(mes.iris, scale=T), main="données normalisées")par(mfrow=c(1,1))Nous utilisons la distance euclidienne
iris.euc <- dist(mes.iris)
iris.scale.euc <- dist(mes.iris.scaled)par(mfrow=c(1,2))
image(t(as.matrix(iris.euc)), main="données brutes")
image(t(as.matrix(iris.scale.euc)), main="données normalisées")par(mfrow=c(1,1))BE et CDA aux points moyensiris.hclust <- hclust(iris.euc)
plot(iris.hclust, hang=-1, cex=0.5)c'est à dire aux options des fonctions dist et hclust
iris.scale.hclust <- hclust(iris.scale.euc)
plot(iris.scale.hclust, hang=-1, cex=0.5)par(mfrow=c(1,2))
plot(iris.hclust, hang=-1, cex=0.5)
plot(iris.scale.hclust, hang=-1, cex=0.5)par(mfrow=c(1,1))iris.scale.max <- dist(mes.iris.scaled, method="max")
iris.scale.hclust.max <- hclust(iris.scale.max)
par(mfrow=c(1,2))
plot(iris.scale.hclust, hang=-1, cex=0.5)
plot(iris.scale.hclust.max, hang=-1, cex=0.5)par(mfrow=c(1,1))iris.scale.hclust.ward <- hclust(iris.scale.euc, method="ward.D2")
par(mfrow=c(1,2))
plot(iris.scale.hclust, hang=-1, cex=0.5)
plot(iris.scale.hclust.ward, hang=-1, cex=0.5)par(mfrow=c(1,1))Les individus dans le plan
iris.scale.kmeans5 <- kmeans(mes.iris.scaled, center=5)
iris.scale.kmeans5K-means clustering with 5 clusters of sizes 12, 50, 37, 24, 27
Cluster means:
Sepal.Length Sepal.Width Petal.Length Petal.Width
1 1.6316667 0.06766667 2.5420000 0.85066667
2 -0.8373333 0.37066667 -2.2960000 -0.95333333
3 0.3863964 -0.20598198 1.0095676 0.37363964
4 0.6858333 0.00100000 1.7503333 0.96316667
5 -0.3137037 -0.43511111 0.1827407 0.01918519
Clustering vector:
[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 5 3 5 3 5 3 5 5 5 5 3 5 3 3 5 3 5 3 5 3 3 3 3 3 3 3 5 5 5 5 3 5 3 3 3 5 5 5 3 5 5 5 5 5 3 5 5 4 3 1 4 4 1 5 1 4 1 4 4 4 3 4 4 4 1 1 3 4 3 1 3 4 1 3 3 4 1 1 1 4 3 3 1 4 4 3 4 4 4 3 4 4 4 3
[148] 4 4 3
Within cluster sum of squares by cluster:
[1] 4.655000 15.151000 11.963784 5.462500 9.228889
(between_SS / total_SS = 93.2 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss" "betweenss" "size" "iter" "ifault"
iris.scale.kmeans5$cluster [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 5 3 5 3 5 3 5 5 5 5 3 5 3 3 5 3 5 3 5 3 3 3 3 3 3 3 5 5 5 5 3 5 3 3 3 5 5 5 3 5 5 5 5 5 3 5 5 4 3 1 4 4 1 5 1 4 1 4 4 4 3 4 4 4 1 1 3 4 3 1 3 4 1 3 3 4 1 1 1 4 3 3 1 4 4 3 4 4 4 3 4 4 4 3
[148] 4 4 3
table(iris.scale.kmeans5$cluster)
1 2 3 4 5
12 50 37 24 27
plot(iris.scaled.acp, col.ind=iris.scale.kmeans5$cluster, choix="ind")Quand une partition est-elle bonne ?
La coupure de l’arbre à un niveau donné construit une partition. la coupure doit se faire :
plot(iris.scale.hclust.ward, hang=-1, cex=0.5)I.intra = numeric(length=10)
I.intra[1] = kmeans(mes.iris.scaled, centers=2)$totss
for (i in 2:10) {
kmi <- kmeans(mes.iris.scaled, centers=i)
I.intra[i] <- kmi$tot.withinss
}plot(1:10, I.intra, type="l")iris.scale.kmeans3 <- kmeans(mes.iris.scaled, center=3)
plot(iris.scaled.acp, col.ind=iris.scale.kmeans3$cluster, choix="ind")heatmap(mes.iris.scaled)my_group=as.numeric(as.factor(substr(variete, 1 , 2)))
my_col=brewer.pal(3, "Set1")[my_group]
heatmap(mes.iris.scaled, RowSideColors=my_col)Mesure de similarité entre deux clustering
à partir du nombre de fois que les classifications sont d'accord
\[R=\frac{m+s}{t}\]
\[ ARI=\frac{RI-Expected RI}{Max RI -Expected RI}\]
cluster.kmeans3 <- iris.scale.kmeans3$cluster
cluster.hclust5 <- cutree(iris.scale.hclust.ward, k=5)
table(cluster.hclust5, cluster.kmeans3) cluster.kmeans3
cluster.hclust5 1 2 3
1 50 0 0
2 0 2 36
3 0 0 26
4 0 24 0
5 0 12 0
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=cluster.kmeans3, choix="ind", title="kmeans en 3 groupes", cex=0.6)
plot(iris.scaled.acp, col.ind=cluster.hclust5, choix="ind", title="hclust en 5 groupes", cex=0.6)par(mfrow=c(1,1))variete <- iris[,5]
table(variete)variete
setosa versicolor virginica
50 50 50
plot(iris.scaled.acp, col.ind=variete, choix="ind")conf.kmeans <- table(variete, cluster.kmeans3)
kable(conf.kmeans)| 1 | 2 | 3 | |
|---|---|---|---|
| setosa | 50 | 0 | 0 |
| versicolor | 0 | 2 | 48 |
| virginica | 0 | 36 | 14 |
variete2 <- rep("notSetosa", 150)
variete2[variete=="setosa"] <- "setosa"
variete2 = factor(variete2)
table(variete2)variete2
notSetosa setosa
100 50
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=variete2, title="variétés observés")
cluster.kmeans2 <- kmeans(mes.iris.scaled, center=2)$cluster
plot(iris.scaled.acp, col.ind=cluster.kmeans2, title="kmeans en 2 groupes")par(mfrow=c(1,1))conf.kmeans <- table(variete2, cluster.kmeans2)
kable(conf.kmeans)| 1 | 2 | |
|---|---|---|
| notSetosa | 97 | 3 |
| setosa | 0 | 50 |
TP <- conf.kmeans[1,1]
FP <- conf.kmeans[1,2]
FN <- conf.kmeans[2,1]
TN <- conf.kmeans[2,2]
P <- TP + FN # nb positif dans la réalité
N <- TN + FP # nb négatif dans la réalité
FPrate <- FP / N # = false alarm rate
Spe <- TN / N # = spécificité
Sens <- recall <- TPrate <- TP / P # = hit rate ou recall ou sensibilité ou coverage
PPV <- precision <- TP / (TP + FP)
accuracy <- (TP + TN) / (P + N)
F.measure <- 2 / (1/precision + 1/recall)
performance <- c(FPrate, TPrate, precision, recall, accuracy, F.measure, Spe, PPV)
names(performance) <- c("FPrate", "TPrate", "precision", "recall", "accuracy", "F.measure", "Spe", "PPV")kable(performance, digits=3)| x | |
|---|---|
| FPrate | 0.057 |
| TPrate | 1.000 |
| precision | 0.970 |
| recall | 1.000 |
| accuracy | 0.980 |
| F.measure | 0.985 |
| Spe | 0.943 |
| PPV | 0.970 |
clues::adjustedRand(as.numeric(variete2), cluster.kmeans2) Rand HA MA FM Jaccard
0.9605369 0.9204051 0.9208432 0.9639434 0.9302767
variete2 <- rep("notVersicolor", 150)
variete2[variete=="versicolor"] <- "versicolor"
variete2 = factor(variete2)
table(variete2)variete2
notVersicolor versicolor
100 50
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=variete2)
cluster.kmeans2 <- kmeans(mes.iris.scaled, center=2)$cluster
plot(iris.scaled.acp, col.ind=cluster.kmeans2)par(mfrow=c(1,1))conf.kmeans <- table(variete2, cluster.kmeans2)
kable(conf.kmeans)| 1 | 2 | |
|---|---|---|
| notVersicolor | 50 | 50 |
| versicolor | 47 | 3 |
TP <- conf.kmeans[1,1]
FP <- conf.kmeans[1,2]
FN <- conf.kmeans[2,1]
TN <- conf.kmeans[2,2]
P <- TP + FN # nb positif dans la réalité
N <- TN + FP # nb négatif dans la réalité
FPrate <- FP / N # = false alarm rate
Spe <- TN / N # = spécificité
Sens <- recall <- TPrate <- TP / P # = hit rate ou recall ou sensibilité ou coverage
PPV <- precision <- TP / (TP + FP)
accuracy <- (TP + TN) / (P + N)
F.measure <- 2 / (1/precision + 1/recall)
performance <- c(FPrate, TPrate, precision, recall, accuracy, F.measure, Spe, PPV)
names(performance) <- c("FPrate", "TPrate", "precision", "recall", "accuracy", "F.measure", "Spe", "PPV")kable(performance, digits=3)| x | |
|---|---|
| FPrate | 0.943 |
| TPrate | 0.515 |
| precision | 0.500 |
| recall | 0.515 |
| accuracy | 0.353 |
| F.measure | 0.508 |
| Spe | 0.057 |
| PPV | 0.500 |
clues::adjustedRand(as.numeric(variete2), cluster.kmeans2) Rand HA MA FM Jaccard
0.53995526 0.07211421 0.07722223 0.57895580 0.40737752
Contact: anne.badel@univ-paris-diderot.fr